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Density functional theories such as the Poniewierski-Stecki theory relate the elastic properties 
of nematic liquid crystals with their local liquid structure, i.e., with the direct correlation function 
(DCF) of the particles. We propose a way to determine the DCF in the nematic state from sim- 
ulations without any approximations, taking into account the dependence of pair correlations on 
the orientation of the director explicitly. Using this scheme, we evaluate the Frank elastic constants 
K\i, K22 and K33 in a system of soft ellipsoids. The values are in good agreement with those 
obtained directly from an analysis of order fluctuations. Our method thus establishes a reliable way 
to calculate elastic constants from pair distributions in computer simulations. 

FAGS numbers: 61.20.Ja, 61.30.Gz, 83.10.Rs 
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I. INTRODUCTION 



Nematic liquid crystals are fluids of anisotropic 
particleS[-,jivhich are aligned preferentially along one 
dircctioru'D. Their orientation is characterized by a di- 
rector n of unit length, with physically identical states n 
and — n. Since the long range orientational order breaks 
a continuous symmetry, the isotropy of space, there exist 
soft fluctuation modes — spatial variations of the director 
n(r) — which cost no energy in the inflnite wavelength 
limit (i.e., the limit where n is rotated uniformlji^ and 
are otherwise penalized by elastic restoring forcesErQ. For 
symmetry reasons, the latter depend on only4;hree mate- 
rial parameters at large finite wavelengthsElt ,, They are 
described by an elastic free energy functional^ 

T{n{r)} = ^ J dr[Kn[y ■ nf + 

K22[n ■ (V X n)]2 + K^^ln x (V x n)]^}, (1) 

which has three contributions: the splay, twist and bend 
modes. The parameters Kaa (a — 1,2,3), called Frank 
elastic constants, control almost exclusively the structure 
and the properties of nematic liquid crystals at meso- 
scopic length scales. Expressions that relate them to the 
microscopic properties of liquid crystals are thus clearly 
of interest. 

Several microscopic agppaaches have been proposed, 
and employed in the pastlZl~E3. Poniewierski^nd SteckiO 
have used the density functional formalismEa to derive a 
set of equations which connects the elastic constants with 
the direct pair correlation function (DCEA- one of the 
central quantities in liquid state theoriestHHEl. In a coor- 
dinate frame where the z-axis points along the director 
n, the equations read 



Ku 



knT 



K22^ 



'xp^^^'{ui^)p^^^'{u2z)uix U2x dv dUi dU2, 
keT 



■ I rl c(r, Ui,U2) 



K: 



'xp^^^'iuiz)p''^^'{u2z)uiy U2y dr dui du2, 
keT 



33 



/ r^c(r,Ui,U2) 



'xp^^^'{uiz)p^^^'{u2z)uix U2x dv dui ciU2, 



(2) 
(3) 
(4) 



where the vector r connects the centers of mass of two 
molecules 1 and 2, ui, U2 are unit vectors along the 
molecule axes, c(r, Ui,U2) denotes the DCF in the ne- 
matic liquid, and p''^^'{uz) is the derivative of the one- 
particle distribution function with respect to Uz. The 
integrals /dr run over all space, and /du over the full 
solid angle, T is the temperature, and fc^ the Boltzmann 
constant. 
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c(r,Ui,U2) 



iS of the form (0)^(0) hatf^ lajter been red- 
and applied in theoriescZl^cZl and simula- 
te study elastic constants in nematic liquid 
The main difficulty with the Poniewierski- 
Stecki equations is that they depend on the DGF in the 
nematic phase, which is not known. Theories have re- 
sorted to approximations, e.g., thajt^use a DCF from an 
effectively isotropic reference sta±££ldl3, or from a state 
with perfectly aligned particlesc3E30. Simulation stud- 
iesHtJ have neglected the explicit angular dependence 
of the pair correlation functions on the orientation of the 
director. Longa et al. have recently pointed out that this 
approximation may not be adequate in nematic liquid 
crystaliiO. 

Alternatively, the elastic constants can also be deter- 
mined directly from the long-wavelength fluctuations of 
the order tensor in Fourier space 



V ^ s 

1=1 



Ui - -I) exp(ik • Vi) 



(5) 



where the sum runs over all particles i in the system, I 
denotes the unit matrix and ® the dyadic product of two 
vectors. The largest eigenvalue of the 3x3 matrix Q 
at zero wavevector (Q(k)|k=o) is the nematic order pa- 
rameter VP2^ and the corresponding eigenvector is the 
director n of the nematic liquid. 

In a reference frame where the z-axis points along n 
and the y-axis is perpendicular to k, the order tensos 
fluctuations have the limiting long-wavelength behaviorQ 



(IQ..(k)| 



2\ fc^o 9 



{P2yVkBT 



4 K^ikl+K^^kl 

k^O 9 (P2fVkBT 



4 K22kl 



K33kl ■ 



(6) 

(7) 



Provided the simulated systems are sufficiently large, the 
elastiie-ceastants can be extracted directly from Eqns. (0), 
(0)0™. 

Allen et al.LiI have used this method to study elastic 
constants in a model liquid crystal, whichj-hiad already 
been investigated earlier by Stelzer et al.E^ using the 
Poniewierski-Stecki equations (||)^(^. The results dis- 
agreed by an order of magnitude. Since the determina- 
tion of elastic constants via Eqns. (g) and (0) is straight- 
forward, it seems reliable and the values calculated by 
Allen et al. ace. presumably accurate. On the other hand, 
Stelzer et al.Lj use an "unoriented nematic approxima- 
tion" , where pair correlation functions are replaced by 
their average over all orientations of the director. Given 
the importance of the Poniewierski-Stecki equations, a 
clearcut test of the applicability of Eqns. (||)^(Q) in a 
nematic liquid crystal is desirable. To the knowledge 
of the present authors, no one has yet employed the 
Poniewierski-Stecki equations with the exact DCF of a 
nematic state. This is presumably due to the fact that 
no method has been proposed so far which allows one to 
extract the full orientation dependent DCF from com- 
puter simulation data. 

The present work attempts to remedy this situation. 
We propose a way to calculate the DCF without any 
approximations from a spherical harmonic expansion 
of the pair distribution function in a uniaxial nematic 
liquid crystal. The expansion coefficients can be de- 
termined from, computer simulations in a straightfor- 
ward mannerEll. A conveniently reformulated version of 
Eqns. (||)^(@) then allows one to calculate the Frank elas- 
tic constants Kn, K22 and ^"33 from a direct inspection 
of expansion coefhcients of the DCF in Fourier space. We 
apply the method to a model system of soft ellipsoidal 
particles in the nematic phase. For comparison, we also 
compute the Frank elastic constants from the fluctua- 
tions of the order tensor, Eqns. (|^) and (0). Wc find that 
the values are in good agreement. Our results thus show 
that the Poniewierski-Stecki theory in combination with 



the correct DCF can be used to bridge between the mi- 
croscopic properties of nematic liquid crystals and their 
mesoscopic, i.e., elastic properties. 

Our paper is organized as follows. We develop the 
theoretical tools needed for our procedure in section II. 
Section III gives details of the simulation model and the 
simulation techniques. The results are presented in sec- 
tion IV and discussed in section V. 



II. THEORETICAL BACKGROUND 

We begin by recalling some common definitionsE3. Let 
us denote by p(u, r) the local number density of parti- 
cles with orientation u at position r. In a uniaxial ne- 
matic liquid at equilibrium with director no, it is dis- 
tributed according to a one-particle distribution function 
(p(u, r)) = p(^^(u), that actually depends on |u-no| only. 
The pair distribution function p(^^(ui, U2,ri — r2) gives 
the probability of finding a particle with the orientation 
Ui at the position ri, and another particle with orien- 
tation U2 at V2- Particles at infinite distance become 
uncorrelated, hence p'-^-'(ui, U2, r) "^Zl^ p(-'^)(ui)p(-'^^(u2). 
This motivates the definition of the so-called total corre- 
lation function 



/i(ui,U2,r) 



p(^)(ui,U2,r) 
p(i)(ui)p(i)(u2) 



1, 



(8) 



which measures the total effect of a particle 1 on a par- 
ticle 2. This effect is often separated into two parts: a 
hypothetical "direct" effect of 1 on 2, characterized by the 
direct correlation function c(ui,U2,r) and an "indirect" 
effect, where 1 is assumed to influence other particles 
3, 4, etc., which in turn affect 2. The total correlation 
functionJs related to the DCF via the Ornstein-Zernike 
equatiocEd 

ft,(ui,U2,ri2) = c(ui,U2,ri2) + 

/ c(ui, U3, ri3) p^'^\vLj) /i(u3, U2, r32)du3dr3, (9) 

where r^ abbreviates r^ — Yj . 

In the framework of density functional theories, the 
direct correlation function has another interpretation as 
the second functional derivative of the excess free en- 
ergy with respect toJocal density distortions 5p{u, r) = 
p(u, r) — p*^^) (u • no)l£j. To lowest order in 5p, the expan- 
sion of the free energy functional about an undistorted 
equilibrium reference state is given by 



5^T- 



knT 



^(ui -U2)(5(ri2) 
p(i)(ui -no) 



-c(ui,U2,ri2; 



X (5p(ui, ri) 6p{vL2,Y2) dri dr2 dui du2. (10) 



In systems of particles with uniaxial symmetry, further 
approximations are not needecEJ. However, the deriva- 
tion is greatly simplified by the additional assumption 



that the relevant long-wavelength distortions can be ex- 
pressed as local distortions of the director n(r), and 
that th£-.density distribution is otherwise at local equi- 
libriumtj 



p(u,r) wp(^^(u-n(r)) 



(11) 



Expanding the free energy in terms of (5n(r) = n(r) — no 
rather than 6p{u, r) and switching to a representation in 
Fourier space, Eqn. (10) then reads 



expressions for the correction terms in systems-ipf asym- 
metric molecules have been given by Yokoyamacj. In this 
paper, we shall be concerned with uniaxially symmetric 
molecules only. 

For practical applications, it is convenient to expand 
all orientation dependent functions in spherical harmon- 
ics Y(„i(u). In the director frame, we obtain 



S^J^- 



VkeT 



(5(ui - U2) 



p(i)(ui -no) 



c(ui,U2,k) 



(12) 



xp^^J'(ui.no)p^^^'(u2-no) 

X [ui • (5n(k)] [u2 • (5n(— k)] rfk dui d\i2. 

This expression has to be related to Eqn. dl^), which has 
the Fourier representation 

.F{n(k)} = i /"dk{i^n[k-n]2 + 

K22[n ■ (k X n)]2 + i^33[n x (k x n)]^}. (13) 



pW(u) = g^/,y,o(u), 

Zoven 

where g is the total bulk number density, and 

F(ui,U2,r)= ^ Fi^jnil2m2lm{r) 



(18) 






Yi,raA^l)Yl2mM2)yim{v). (19) 
F(ui,U2,k)= ^ Fi^rnil2m2lrn{k) 



ll,l2.' 

711, 7712, TT 



To this end, we expand the DCF c(ui, U2, k) in Eqn. (|l^ 
in powers of k up to second order. For convenience, we 
choose a coordinate frame such that the z-axis points in 
the direction of no (director frame). 

Since a global rotation of the director n does not 
change the free energy, the leading term k = must 
vanish, i.e., one has 



rz,™,(ui)ri,™,(u2)i,™(k). (20) 

Here F stands for any of p^'^\ h, or c, f denotes the unit 
vector r/r, and k the unit vector k/fc. The symmetry 
of the nematic phase ensures that all coefficients are real 
and only coefficients with m + mi +7712 = 0, and I + I1 + I2 
even, enter the expansions (|9|) and (|2^). If the molecules 
have uniaxial symmetry, every single li has to be even in 
addition. 

Next we derive matrix versions of Eqns. (|8[) and (|^). 
To simplify the expressions, we introduce the notation 



P^'H^z) 



ldu= /"c(ui,U2,k = 0)p(l)'(7/i,,) 

^^^'{u2,z)ui,aU2.adUidU2 (14) 



r: 



U'l" 



JduY*Ju)Yi,^„A^),Yi^',m"{^) (21) 



X P 



for a = x,y. Eqn. |£j) has been derived in a different 
context by Gubbina_3 and is quite generally valid. For 
symmetry reasons, the terms linear in k in the expan- 
sion of ( |l^ ) vanish too. The quadratic terms lead to an 
expression of the form (O) , with Ka given by 



Kii = - 



ksT /•c>2c(k,Ui,U2) 



^'^^"^^J}'^^'t^^ C{l"l'l; 000)Cil"l'l; m"m'm), 
Att[21 + 1) 



where C are the Clebsch-Gordan coefficients. The total 
correlation function h can then be calculated from p^^' 
by inversion of the matrix version of Eqn. (|^) 



n(2) 



k=0 



K22^- 



2 J dkl 

X-P^^^'{ui^)p^'^^'{u2z)ui^ U2x dVLi dVL2 

keT /•c>2c(k,ui,U2) 



dkl 



k=0 



^33 = - 



Xp^ > {Ui^)p^ > {u2z)uiyU2ydUidU2, 

keT f d\{'k,u,,U2) 



(15) 



(16) 



PhUl2m2lmir'> = ^ [^'iT^fllfl2Sm,oSm2oSloSrnQ 



dkl 



k=0 



'XP^^^'{ui^)p'^^'^'{u2z)uix U2x dUi (iU2, (17) 

which is the Fourier space version of the Poniewierski- 
Stecki equations (||)-(||). As mentioned above, the same 
result can be derived without the approximation (O) for 
systems of particles with uniaxial symmetrycS. Compact 



(22) 



Eqn. (p2) is a linear system of equations and can be 
solved for the coefficients of h by standard numerical 
methods. 

The Ornstein-Zernike equation (O) is most conve- 
niently solved in Fourier space k. We calculate the co- 
efficients hi-^mii2m2im{k) of thc total correlation funC|tipn 
in Fourier space by using thc Hankel transformationO 



hi^mii2m2im{k) = ini' r'^ ji{kr) hi^mii27n2i7n{r) dr, (23) 
Jo 



with the spherical Bessel functions ji. The matrix ver- 
sion of the Ornstein-Zernike equation (j^ in Fourier space 
reads 

'^l-i_niil2m2lm\'^ ) — ^Iini-Ll2m2lm\/^ ) 



l'm'l"m"i. 



^ Jig -■- mm'm"-'- rjJamsOV -■-/ ■ V^^^ 



The result for the direct correlation function c(fc) is read- 
ily transformed back into real space by another Hankel 
transformation. However, this is not necessary for our 
purpose, because the Poniewierski-Stecki equations as- 
sume a very handy form in Fourier space: the spherical 
harmonic representation of Eqns. ([l5[) ~ (|l7|) reads 



^""2dF^"^^^ 



fc=0 



for i^ 1,2,3 



(25) 



with 



Cuik) 



Y. Vhill + l)Vl2{l2 + l) fl, h 



hh 



8v^ 



I [Qii/2-ioo(fc) + Qi-i;2ioo(fc)] 

+ 'Vi — [ci^U2~i2o{k) +c;i-ii2i2o(fc)] 



+ w^—r=-[ci,ii^i2-2{k) + ci,-ii^-i22{k)] I (26) 

and {vi,V2,vz) = (-1,-1,2), {wi,W2,W'i) = (-1,1,0). 
Deriving these equations, we have exploited the relation 
JdrF{r)r^ = — 9^_F(k)/(9fc^|k=o and properties of spher- 
ical harmonics. Finally, Eqn. dl^) can be rewritten as 

Cuik = 0) = -keTTT / du, (1 - ul) \..) '[ , (27) 

J-l p'^'-'iUz) 



where Cii{k) is defined as in Eqn. (|2i 



III. MODEL AND SIMULATION DETAILS 

We performed computer simulations of a system of axi- 
ally symmetric rigid particles, which interact via a simple 
repulsive pair potential 



V^, 



4eo {Xlf - Xg.) + eo : Xf^ > 1/2 
: otherwise 



Here Xij — oa/ijij — Gij + gq), rij denotes the distance 
between particles i and j , and the shape function 



cry(Uj,Uj,fy) =cro i 1- o 

(Ui • fij 



X 



(Uj -Tij +Uj -Tij^ 
1 + XU., • Uj 
-1/2 



1 - XUj • Uj 



(29) 



approximates the contact distance between two ellipsoids 
of elongation k = cre„d-e„d/CT.ide-sido = V(l + x)/(l-x) 
with orientations u^ and u^, which are separated h^-a 
center-center vector in the direction of f^ — Tij/viji^. 
We use throughout scaled units defined in terms of eo, 
(To, the particle mass mo and the Boltzmann constant fc^. 
We studied systems of particles with elongation k = 3 at 
temperature T = 0.5 aud number density g = 0.3. The 
pressure was P — 2.60ll3. This corresponds to a state 
well in the nematic phase: at fixed temperature T = 0.5, 
the fluid remains nematic down to the density g — 0.29 
or, equivalently, the pressure P — 2.35E3. The average 
order parameter density in our system was (P2) = 0.69 
and the fourth rank parameter was {P4{u ■ no)) ~ 0.31, 
Pi{x) = (35a;'* — 30a:^ -I- 3)/8 being the fourth Legendre 
polynomial. 

The pair distribution function was determined in sys- 
tems of TV = 1000, 4000 and 8000 particles in cubic 
boxes with periodic boundary conditions. For the N = 
1000 system we used a Monte Carlo (MC) program by 
H. Langgl3. Trial moves picked a particle at random and 
attempted in random order either a rotation or a trans- 
lation, with maximum step sizes chosen such that the 
Metropolis acceptance rate was roughly 30%. The larger 
systems were studied with a massively parallel computer, 
using a domain decomposition molecular dynamics (MD) 
program, that has been codeveloped by one of us (GG). 
These simulations were performed in th e jni crocanoni- 
cal ensemble usixig the rattle integratoiOlIa with time 
step At = 0.003O and molecular moment of inertia / = 
2.5. Run lengths were 8 million MC steps, one MC step 
consisting of 2A'' trial moves, or 10 million MD steps, re- 
spectively; data for the pair distribution function were 
collected every 1000 or 10000 steps. 

The order tensor fluctuations are sampled most effi- 
ciently if the k-vectors in Eqns. (g) and (0) are always 
on the same grid. They were therefore determined from 
independent simulations in an ensemble where the direc- 
tor Mf) was constrained to the Z-axis of the simulation 
boxO. Thus the xyz frame of Eqns. (g) and (^ becomes 
coincident with the XYZ frame of the simulation box. 
The constraint was implemented in the MD simulations 
by adding two global Lagrange multipliers to the inte- 
grator, so that Qxz{0) — Qyz{0) = at every time 
step. Our procedure was similar to that introduced by 
AlLea et al.O, except that we used an improved integra- 
torEy designed in the spirit of rattleOIIS, so that it is 
symplectic and fulfills the constraiats exactly. The same 
integrator has already been usedEa to calculate K22 in 
a Gay-Berne fluidEd; the value compared well with an 
estimate from a thermodynamic perturbation approach. 
Here, we simulated a system of A^ = 4000 particles in 
a cubic box over 10 million MD steps, and a system of 
N = 16000 particles in aa elongated box with side ratios 
Lx : Ly : Lz = 1:1: 2e3 over 5 million MD steps. Data 
for the order tensor were collected every 200 steps. The 
largest autocorrelation times were of the order of 10^ MD 
steps at the lowest k values and dropped rapidly below 
1000 MD steps for higher k. 



IV. DATA ANALYSIS AND RESULTS 

We begin by presenting the results for the order tensor 
fluctuations. FoUowing Ref. ^ we calculated the quan- 
tities 



W..(k) 



W,.(k) 



9{P2fVkBT k- 
4(|Q..(k)P) 

^{P2?VkBT k- 



Kiikl + Ks^kl 



K22kl + ifsafc^ 



(30) 



(31) 



4(|Q,.(k)P) 

where the frame is chosen such that k lies in the xz- 
plane (cf. Eqns. (||) and (0)). More specifically, we 
evaluated the order tensor in Fourier space Q(k) on 
a k-grid with 6x6x6 grid points in the small sys- 
tem {N = 4000), and 6 x 6 x 12 grid points in the 
large system [N — 16000). Then we applied a rotation 
Q(^y^)(k) = U(k)Q(^'^^)(k)UT(k) into the desired co- 
ordinate frame such that ky ~ 0, and calculated the aver- 
ages (IQQz(k)p) and the Wqz (k)-surface in that frame. 
Because of the constraint on no, U(k) is a constant 
throughout the run. 

In the high wavelength-limit k -^ cx3, Wazi}^) (q:=1,2) 
takes the valueO 



>Va.(k) 



{P2)^pkBT 



((P2)/21-4(F4)/35 + l/15- 



(32) 




FIG. 1. W^, (a) and Wyz (b) surfaces for iV = 4000 (cubic 
box) and A*' — 16000 (elongated box); the smaller, finer spaced 
grids correspond to the larger system. The fits (dotted lines) 
coincide almost perfectly with the data (solid lines). 



In our simulations, we obtained 1.13, which is in good 
agreement with the theoretical value 1.12. 

The results for the Waz (k) surfaces are shown in Fig. [|. 
The data for the small system (coarse grid) match almost 
exactly those for the large system (fine grid). They were 
fitted to a fourth order polynomial in k^ and fc^ (i.e., 
with highest order terms fc^, /c^fc^, • • • , /cf ) without a ze- 
roth order term. Higher orders were disregarded because 
the 4th order coefficients turned out to be already very 
small. Normal equations and singular value decompo- 
sition gave the same results. Fig. H demonstrates that 
the fit is almost perfect. The leading coefficients give the 
elastic constants, shown in table |. As expected for elon- 
gated molecules, one finds that K33 is largest, followed 
by Kii and K22- 

Next we discuss the results for the pair correlation 
functions. The spherical harmonics expansion coeffi- 
cients of theLjpair distribution function p'^' were deter- 
mined usingEil 



rlimil2m2lm\ ) 



Y, 



l\m 



An g^ g{r) 
^{u,)Y*^^{u2)Y,:^{r))sr 



(33) 



where {■)sr denotes the average over all molecules in a 
shell Sr from r to r + Sr, and the function g{r) is the 
number of molecular centers at distance r from a given 
molecular center, divided by the number at the same dis- 
tance in an ideal gas at the same density. The calculation 
of these averages is very time consuming, since a great 
number of coefficients has to be evaluated, and was there- 
fore carried out in part on a massively parallel machine. 
We have determined coefficients for values of IJi up to 
Unax = 6 in all systems, and for values up to Imax — 8 
in the smallest system. The bin size was 6r = 0.04 and 
the cutoff distance r^ax was chosen to be 40,% of the box 
side L in order to reduce boundary effectg23. 

From the pair distribution function we calculated the 
total correlation function by inverting Eqn. (E2). The lat- 
ter was then Fourier transformed according to Eqn. (|23|). 
There is a subtle problem here: due to the elasticity of the 
nematic phase, the total correlation function decays al- 
gebraically like 1/r. This follows directly from Eqn. dfi)Q. 
Before applying Eqn. (Eq), we thus fitted the simulation 
data points at the largest distances r > tq to a power law 
of the form b/r and extrapolated h{r) to infinityE3. The 
parameter ro was chosen to be 2.8, 4.0 and 5.3 in systems 
of TV = 1000, 4000 and 8000 particles, respectively. 



System 
size 



4000 
16000 



Order tensor fluctuations 

(Kii) {K22) {K33) 



0.53 ± 0.01 
0.53 ± 0.01 



0.30 ± 0.01 
0.30 ± 0.01 



1.60 
1.59 



0.01 
0.01 



TABLE I. Elastic constants from the analysis of order 
tensor fluctuations for systems of different size N. 
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-20 
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Imax = 6 and Imax = 8 do not differ significantly. 

Since the calculations with Imax = 8 were very 
time consuming (one has 1447 different expansion coef- 
ficients), we used Imax = 6 in the analyses of the larger 
systems (469 different expansion coefficients). 

The results were the same for systems of size A^ — 
1000, 4000 and 8000. Furthermore, they were not affected 
by the presence of a director constraint: as mentioned in 



10 



15 



Cn(k) 



• ---• N=8000 

<^ N=4000 
* * N=1000 



FIG. 2. Expansion coefficient /i2i2-i2o('') of the total cor- 
relation function h vs. r in systems of size A*' — 8000 (solid 
line), N = 4000 (dotted line), and TV = 1000 (dashed line). 
Cutoff radii were rmax = 11.9,9.4, and 6.6, respectively. The 
long dashed line indicates the extrapolation towards r —> ca 
(for the dataset A'' = 1000). Inset shows same data vs. 1/r. 

It turned out that the long-range tail was quite pro- 
nounced for coefficients of h with mi = ±l,r7i2 = ±1, 
and almost negligible for the others. In Fig. || we show 
an example of a coefficient with a pronounced long-range 
tail, the coefficient hi-^mii2m2im{r) with li = I2 = I = 2, 
mi — 1, m,2 = — 1 and m, = 0. The data for different 
system sizes N = 1000, N = 4000, and A^ = 8000 he 
almost on top of each other, hence the form of h{r) at 
T < fmax is not affected by noticeable finite size effects. 
The dominating finite size problem comes from the un- 
certainty of the extrapolation, if the available range of 
h{r) is too short. 

The rest of the analysis was straightforward. From the 
coefficients of the total correlation function in Fourier 
space, hi^mii2m2im{k), those of the DCF were obtained 
by solving the linear matrix equation (|2j) . Then we cal- 
culated the functions Cii{k) as defined in Eqn. (Eq). Ac- 
cording to Eqn. (psj), the elastic constants Ku can be 
determined from the initial slopes in a plot of Cu(k) ver- 
sus k^ . Data for Cii{k) are shown for different system 
sizes in Fig. (P). The points at zero wavevector Cm(0) 
were calculated using Eqn. (p7|). They fit nicely on the 
straight lines at fc -^ 0, hence the data are consistent 
with the requirement (O) or (CT)- This gave additional 
confidence in the quality of the analysis. The slopes of 
the straight lines yield the elastic constants. 

The results are summarized in table ||. We have cal- 
culated the DCF from the pair distribution function p^^' 
using an upper cutoff Imax = 2,4 and 6, respectively, in 
the matrix equations (p^) and (p4|). Already the lowest 
order calculation with Imax — 2 gave elastic constants of 
the correct order of magnitude. Quantitatively reliable 

resuhs were obtained with Imax > 6: we checked in the See section |][V| for details, 
smallest system that the results from calculations with 
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FIG. 3. Weighted sum of the DCF expansion coefficients 
Cii{k) as defined in Eqn. (M) vs. k^ for different system sizes 
A^ (unconstrained director, evaluated using coefficients up to 
Imax = 6). The points at fc = are taken from Eqn. (pTI). 
The initial slopes give the elastic constants Ku. Thick solid 
lines indicate corresponding fits for the A' = 4000 system. 



System 




Direct correlation function | 


size 


Z„,ax 


{Kii) 


{K22) 


{K,s) 


1000 


8 


0.55 ± 0.02 


0.35 ± 0.03 


1.56 ± 0.04 




6 


0.51 ± 0.02 


0.34 ± 0.03 


1.52 ± 0.04 




4 


0.53 ± 0.03 


0.23 ± 0.02 


1.32 ± 0.04 




2 


0.51 ± 0.01 


0.20 ± 0.01 


1.56 ± 0.04 


4000 


6 


0.51 ± 0.02 


0.31 ± 0.01 


1.51 ± 0.03 




4 


0.65 ± 0.02 


0.27 ± 0.02 


1.23 ± 0.03 




2 


0.53 ± 0.01 


0.22 ± 0.01 


1.46 ± 0.03 


*4000 


6 


0.52 ± 0.02 


0.31 ± 0.01 


1.51 ± 0.03 




4 


0.65 ± 0.02 


0.27 ± 0.02 


1.24 ± 0.04 




2 


0.53 ± 0.01 


0.22 ± 0.01 


1.48 ± 0.03 


8000 


6 


0.51 ± 0.02 


0.33 ± 0.02 


1.48 ± 0.03 




4 


0.61 ± 0.01 


0.29 ± 0.02 


1.25 ± 0.04 




2 


0.54 ± 0.01 


0.23 ± 0.01 


1.47 ± 0.04 



TABLE II. Elastic constants from the DCF method for 
systems of different size A'. (*) marks a system that has 
been simulated with a director constraint. Results are shown 
for different choices of the cutoff value Imax in the spherical 



harmonics expansion of the pair distribution function p 



(2) 



section [II, the DCF was mostly calculated in uncon- 



strained systems, but we also studied the DCF in one 
constrained system for comparison. 

Finally, we compare the values of the elastic constants 
calculated by the DCF approach with those obtained 
from the order fluctuation analysis, shown in table y. The 
values for Ku and K22 are identical for both methods. 
K33 is slightly underestimated by the DCF analysis with 
Imax — 6, but the result increases with Imax, and agrees 
within the error with that of the order fluctuation anal- 
ysis at Ijnax = 8. 

One might ask how much the successive coefficients of 
the (correct) DCF contribute to the elastic constants. We 
found that the contribution of the coefhcients with /i > 4, 
I2 > 4: or I > 4 is very small. If we include only terms up 
to /,/j = 4 in Eqn. (p6|), we obtain Ku = 0.55, K22 = 
0.21 and ^^33 = 1.51 (in the largest system N = 8000), 
which is very close to the final values quoted in table ||. 
However, we could not push this analysis further. If we 
include only terms up to I, k — 2, the resulting Cu{k) are 
very concave and have no well-defined initial slope in a 
plot vs. fc^. Hence the contributions of successive Is to 
the elastic constants cannot be distinguished. 



V. SUMMARY AND CONCLUSIONS 

We have presented a method which allows one to deter- 
mine without approximations the direct correlation func- 
tions in nematic liquid crystals from computer simula- 
tions, and to calculate elastic constants op-that basis ac- 
cording to the Poniewierski-Stecki theoryEJ (||)^(^). We 
have applied this method to a nematic fluid of soft el- 
lipsoids. In the same system, the elastic constants were 
also determined by an established approach, the analysis 
of order tensor fluctuations. 

Our study represents a direct test of the Poniewierski- 
Stecki theory. We found that the results obtained with 
the two methods agree well with each other. The 
Poniewierski-Stecki theory can thus be employed to cal- 
culate elastic constants, at least in our system, provided 
that the exact direct correlation functions are used in the 
equations. 

Hence we have established an alternative way of calcu- 
lating elastic constants in nematic liquid crystals. As 
long as a simulation is performed solely to determine 
elastic constants, the order tensor fluctuation approach is 
still more efficient: the statistical error of pair correlation 
functions must be quite small for a reliable DCF analysis, 
and the analysis is very time consuming. However, the 
DCF approach has the advantage of being straightfor- 
ward; elastic constants can be computed from arbitrary 
bulk simulations, if the pair distribution functions are 
known with sufhcient accuracy. Even the calculation of 
spatially varying elastic constants, e.g., in the vicinity of 
surfaces, is conceivable. 

The direct correclation function is a central quantity 



in liquid state theories. The study of direct correlation 
functions in the nematic phase is therefore interesting 
in its own right. We shall examine them in more detail 
and compare them to tiose in the isotropic phase in a 
forthcoming publicationEJ. 
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